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A few years ago, Bouchaud et al. introduced a phenomenological model to describe surface 
flows of granular materials [J. Phys. Fr. I 4, 1383 (1994)]. According to this model, one can dis- 
tinguish between a static phase and a rolling phase that are able to exchange grains through an 
erosion/ accretion mechanism. Boutreux et al. [Phys. Rev. E 58, 4692 (1998)] proposed a modifi- 
cation of the exchange term in order to describe thicker flows where saturation effects are present. 
However, these approaches assumed that the downhill convection velocity of the grains is constant 
inside the rolling phase, a hypothesis that is not verified experimentally. In this article, we therefore 
modify the above models by introducing a velocity profile in the flow, and study the physical con- 
sequences of this modification in the simple situation of an avalanche in an open cell. We present a 
complete analytical description of the avalanche in the case of a linear velocity profile, and generalize 
the results for a power-law dependency. We show, in particular, that the amplitude of the avalanche 
is strongly affected by the velocity profile. 
Short title: Thick granular surface flows 
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I. GENERAL PRINCIPLES 



A. Onset of avalanches 

It is a dailylife experience that the top surface of a mass 
of granular matter need not be horizontal unlike that of a 
stagnant liquid. However, there exists an upper limit to 
the slope of the top surface, and the angle between this 
maximum slope and the horizontal is known, for non- 
cohesive material, as the Coulomb critical angle 9 max . 
Above this angle, the material becomes unstable, and an 
avalanche at the surface might occur. The Coulomb angle 
is related to the friction properties through tan 9 max = fa 
where fii is an internal friction coefficient [ [l]] . 

As of today, the physical picture associated with the 
onset of the avalanche is still obscure. One could imagine 
a local scenario in which the dislodgement of some un- 
stable grains leads by amplification to a global avalanche 
(see for instance [§]). Alternately, one can think of a 
delocalized mechanism [ |j| , in which a thin slice of mate- 
rial is destabilized and starts to slide as a whole. In the 
present paper, we will focus on the latter point of view. 

It has been recently suggested [ || that the thickness 
of the initial gliding layer should be of the order of £, 
the mesh size of the contact force network [ ||-|(| . For 
simple grain shapes (spheroidal), one expects £ ~ 5-10 
grain diameter d. The angle at which the avalanche pro- 
cess actually starts is of the order of 9 max +£/L, where 
L is the size of the free surface. At the moment of onset, 
our picture is that this initial layer starts to slip, and 
is rapidly fluidized by the collisions with the underlying 
heap, therefore generating a layer of rolling grains on the 



whole surface. 

Now that we have proposed a description of the initial 
situation, we may turn to the model scheme accounting 
for the further evolution of the avalanche. 



B. Saturation effects for thick avalanches 

Some years ago, Bouchaud, Cates, Ravi Prakash and 
Edwards introduced a model to describe surface flows 
of granular materials [ 0] . The model assumes a rather 
sharp distinction between immobile particles and rolling 
particles and, accordingly, introduces the following two 
important physical quantities (see Fig. |l|): the local 
height of immobile particles h(x, t) (where x denotes the 
horizontal coordinate [ |) and t the time), and the local 
amount of rolling particles R(x, t). 

The time evolution of h(x, t) is written in the form 



dh 

at 



= lR(0n ~ 



(1) 



where 9 ~ tan(0) = dh/dx is the local slope, 7 a char- 
acteristic frequency and 9 n the neutral angle of grains at 
which erosion of the immobile grains balances accretion 
of the rolling grains. For the rolling particles, Bouchaud 
and co-workers wrote a convection-diffusion equation [ [?]] 
that was later simplified by de Gennes as [ ||] 



dR _ dR dh 
dt dx dt 



(2) 



where v is the downhill typical velocity of the flow, and 
is assumed to be constant. 
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According to the Bouchaud-Cates-Ravi Prakash- 
Edwards (BCRE) model, dh/dt is linear in R [see 
Eq- (0)]- This is natural at small R, when all the rolling 
grains interact with the immobile particles. But as ex- 
plained in Refs. [E0,O], this cannot hold when R be- 
comes larger than a given saturation length since the 
grains in the upper part of the rolling phase are no longer 
in contact with the immobile grains. The length £' is 
expected to be of the order of a few grain diameters 
d [|l^]. This led Boutreux, Raphael, and de Gennes to 
propose [ [TT| a modified saturated version of the BCRE 
Eq. (Pj), valid for thick surface flows and of the form 



dh to 



(r » a 



(3) 



where v u h is defined by v u h = 7^'- The constant v u h has 
the dimensions of a velocity. 

The description of thick avalanches modelized by 
Eq. (||) was discussed in Ref. [ pH|. However, one might 
encounter situations where the local amount R of rolling 
particles is rather large except in some regions of space 
where it takes values smaller than For such cases, 
various 'generalized' forms of the BCRE equations valid 
both in the large and small R limit, and able to han- 
dle intermediate values have been proposed [ [l(i| . [l3| . [l^| . 
As we will be concerned only with thick flows, we will 
henceforth use the saturated form (H). 



C. Velocity profiles in thick flows 

We now consider the hypothesis made in Eq. (0) that 
the downhill typical convection velocity of the rolling 
grains v is constant. As a matter of fact, v might vary 
for two reasons. 

First, v depends on the local slope dh/dx of the static 
bed, reflecting that the mean convection velocity should 
increase as the sandpile is further tilted. However, in 
the situation we are going to consider, the slope should 
never depart from 8 n by more than a few degrees, so that 
the variations of v originating in this may reasonably be 
taken to be negligible. 

Second, v might as well depend on the local amount 
of rolling particles R. This dependence is quite natu- 
ral, since as soon as the thickness of the flow exceeds a 
few grain diameters, one would expect a velocity gradient 
perpendicular to the flow to establish. Such a possibility 
was already considered by Bouchaud et al. [ [L4|], but, to 
our best knowledge, not further studied. We think that 
taking this velocity gradient into account does lead to 
an improvement of the model description of avalanches. 
In the forthcoming sections we will analyse the physical 
consequences of this modification. 

If analyticity is assumed, we can expand v(R) in pow- 
ers of R, and considering only the first two contributions 
to be significant, we write: 



with r a constant, homogeneous to a shear rate, and vq 
a constant velocity. 

When R becomes small, Eq. (Q) tells us that v(R) be- 
comes constant [v(R) — > vo}. Physically this velocity 
should correspond to the typical convection velocity of 
a single grain on a bed of immobile grains. For simple 
grain shapes (spheroidal) and average levels of inelastic 
collisions, one expects this velocity vq to scale as (gd) 1 / 2 
(where g is the gravity) [ pj . Similarly, the shear rate T is 
expected to scale as (gd) ^/d ~ (g/d) 1 ^ 2 [ We can 
therefore rewrite Eq. (Eh: 



v(R) = T(R + d). 



(5) 



We note that Vq becomes negligible compared to TR as 
soon as R exceeds a few grain diameters. 

In our approach, the typical velocity v(R) depends 
linearly on the local rolling height R [Eq. (g)]. Such a 
form is in part motivated by the r ecent work of Douady 
et al. [ |l6| (see also Section IV C ) . It is also supported 
by the experimental results of Rajchenbach et al. who 
carried measurements in a rotating drum [ |l7],D . These 
authors have found linear velocity profiles in the surface 
flow, with a shear rate T independent of the thickness of 
the flow. However, in other experiments of chute flows 
carried on rough inclined planes, Anzaza et al. [ [l8| and 
Pouliquen [ observe that the mean velocity (averaged 
on cross-sections) scales as a power-law of the thickness 
with an exponent about 3/2. In the following we will 
mainly focus on the linear form (EJ), since it allows us 
to give explicit analytical solutions, and shall discuss the 
changes that are to be b rought in the case of a power-law 
velocity in Section IV A. 

In the next section, we will derive the governing equa- 
tions from the saturated BCRE equations and the above 
considerations on the velocity profile inside the flow. 



D. Governing equations 

We may define a reduced profile h, deduced from h by 
substracting the 'neutral' profile 9 n x: 

h(x,t) = h(x,t) — 9 n x. (6) 

Using Eqs (||), (^|), (j|) and (|6|), we easily obtain the 
following system: 



dh 
~di 



Vuh 



dh 
d.r 



OR ^dR dh 



(7) 
(8) 



v(R) =v Q + TR. 



(4) 



In our approach, Eqs. (0) and (||) are the governing 
equations for surface avalanches displaying linear velocity 
profiles. 

An important point is that we must have R > for 
Eqs. (0) and (|) to be valid. If we reach R = in a 
certain spatial domain, then Eq. (0) must be replaced in 
that domain by dh/dt = 0. 
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II. APPLICATION TO THE SIMPLE CASE OF 
AN OPEN SYSTEM 

A. Physical situation 

We will now solve Eqs. (]?]) and (||) in the following 
simple situation: we consider a cell, of dimension L, par- 
tially filled with monodisperse grains of diameter d, as 
shown on Fig. |[ The heap has an initial uniform slope 
9max, the Coulomb angle of the material. The origin of 
the x axis is taken at the bottom of the cell, and the 
orientation of the axis is such that the slope of the heap 
is positive. 

We now consider that an avalanche has just started in 
the cell (see Section I A ) , so that we have at time t = a 
layer of rolling grains in the whole cell, of thickness ~ £ 
greater than the saturation length We may thus use 
the saturated equations (Q) and (||) from the beginning 
of the avalanche. 

As the rolling population will rapidly grow and become 
independent of the initial thickness £ (for £ small) , we can 
as well consider the initial condition on R to be: 



(where the subscript uh stands for 'uphill'). 

The wave starts from the bottom of the cell at time 
t = and reaches the upper end at time t 2 defined by: 



R{x,t = 0) = 0. 
We also know the initial value of h: 

h(x, t = 0) = (9 max -6 n )x = r]X, 



(9) 



(10) 



where 77 is defined as the (positive) difference between 
the Coulomb angle and the neutral angle. 

We have additional conditions in our system, due to 
the boundaries. At the top of the cell, there is no feeding 
in rolling species, so that we impose: 

R(x = L, t) = at any time t > 0. (11) 

Another condition arises from the fact that grains fall off 
the cell at the bottom and cannot accumulate there: 

h(x = 0, t) = at any time t > 0. (12) 



B. Uphill wave in the static phase 

Equation ^ can be readily solved along with condi- 
tions (0) and (H) to give: 

h(x, t) = rjv u h,H{x — v u ht) x for < x < L, (13) 

where H denotes the Heavyside unit step function 
[H{u) = 1 if u > 0, H(u) = otherwise]. This re- 
sult corresponds to the uphill propagation (at constant 
speed v u h) of a surface wave on the static phase. Let us 
call x u h{t) the time-dependent position of the wavefront, 
given by: 



t 2 = L/v uh . 



(15) 



x u h(t) = V u ht. 



(14) 



At a given time t (smaller than £2), the profile of the 
static phase can be described as follows: ahead of the 
wavefront [x u h{t) < x < L], the profile is linear and the 
slope is the initial angle 6 max (since h = r]v u h)- Behind 
the wavefront [0 < x < x u h(t)], the slope has decreased 
and reached the neutral angle 9 n (h — 0) (see Fig. ||). 
For times t > t% the slope of the static phase inside the 
cell is uniformly equal to the final value 8 n , which is thus 
the angle of repose of our specific open cell system [ po) . 



C. Downhill convection of rolling grains 

Substituting Eq. (|l^) into the evolution equation (||) 
for R gives: 

dR dR 

— -T{R + d)—=i 1 v uh H{x-v uh t) (16) 

Eq. ( |l6| ) is a non-linear convection equation. The 
rolling species are thus convected downhill, with a con- 
vection velocity dependent on the local rolling thickness 
R. In the spatial region x > v u ht 1 the right-hand side 
(which couples the evolution of R to that of h) plays the 
role of a source term, leading to an amplification of the 
avalanche. On the contrary, for x < v u ht, the right-hand 
side goes to zero, so that the material flowing through the 
surface x = v u ht from uphill is simply convected, without 
amplification nor damping. 

Equation (|l6| ) can be solved analytically by using the 
method of characteristics [ |2l| , which utilizes the prop- 
erty that certain types of partial differential equations 
reduce to a set of ordinary differential equations along 
particular lines, known as the characteristic curves. For 
more details on this method and on its application in the 
case of Eq. (p^), see Appendix. 



D. Propagation of boundary effects in the cell 

Before we go to the precise solutions, we can try to 
get some physical insight of the way the avalanche is go- 
ing to develop. The global shape of the rolling phase 
at different moments during the avalanche is of course 
very dependent on the boundary condition (|ll]) for R in 
the cell, but also on the condition (|l2|) for h, since the 
evolution of h and R are coupled. 

However, the effects of these boundary conditions can- 
not spread over the entire cell instantly after the be- 
ginning of the avalanche, and shall propagate with fi- 
nite velocities. We then expect the progression of these 
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boundary effects (one could say, the propagation of the 
'information' on the boundaries) to control the evolution 
of both the rolling and the static phase. For instance, 
in the case of the static profile h, Eq. (|l3| ) tells us that 
the bottom boundary condition (Jl^) [h(x = 0, t) = at 
any time t] brings progressively h to zero everywhere in 
the cell, and also, that the propagation proceeds with a 
velocity v uh - 

Hence, we expect that the description of the avalanche 
should naturally split up in different temporal 'stages', 
according to the degree of extension of the different 
boundary effects, and that the cell should divide in sev- 
eral 'regions', according to whether it is under the in- 
fluence of the top boundary condition or the bottom, or 
both, etc. This shall become clear as we will now go into 
the precise description of the avalanche. 



III. UNFOLDING OF THE AVALANCHE 



A. Stage I: The avalanche grows to maturity 

This stage starts at t = with the beginning of the 
avalanche. From the above considerations, we know that 
the boundary effects start to propagate with finite veloc- 
ities from both ends of the cell. We can therefore define 
'propagation fronts' for these effects: we call Xdh(t) the 
position of the front originating in the boundary con- 
dition at the top of the cell [the subscript dh means 
that the motion of this front is downhill], and x u h(t) 
the corresponding 'uphill' front, originating in the bottom 
boundary condition, and that we already defined earlier 
(t) = v u ht [Eq. (Q)]. Figure ^-a presents a typical 
picture of the situation during Stage I, where the fronts, 
after leaving their respective cell ends, move in opposite 
directions and one toward the other. As a consequence, 
they shall finally meet at a certain time, that we hereafter 
denote t\. This time t\ defines the end of what we call 
'Stage F (which is thus characterized as the time interval 
< t < t%), and the beginning of 'Stage IF (described in 
next section). 

The relative positions of the fronts naturally define 
three spatial regions in the cell (Fig. f|-a). To the left 
of x u h, the effects of the bottom boundary condition 
[Eq. (|l2|)] are predominant. We call this region the bot- 
tom region. We remark that it is constantly extending 
uphill during Stage I [following the motion of x u h(t)]. To 
the right of Xdh(t), we define the top region, which ex- 
tends downhill, and where the evolution of the avalanche 
is controlled by the upper end condition [Eq. (jll|)]. Fi- 
nally, between those two regions remains a central region, 
where none of the boundary effects can yet be felt. This 
last region shrinks during Stage I, and ultimately disap- 
pears at time t — t\ when the bottom and top region 
connect. We now describe the precise evolution of the 
avalanche, region by region [we will only give the form of 



the rolling amount R(x,t), since h(x,t) is already known 
from Eq. ©]. 



1. Top region 

From the above definition, the top region corresponds 
to the spatial domain Xdh(t) < x < L. Within this do- 
main, Eq. (n6f) reads: 



m- T{R + d) dx- 



V v uh 



(17) 



[since x > x dh (t) > x uh = v uh t]. 

Solving this equation with the boundary condition 
R(x — L, t) — (see Appendix for details) gives the ex- 
pression of R valid in this region: 



R(x,t) = -d+ J<P + 2(L - x 



V V uh 



(18) 



We also obtain the precise position of the 'downhill front' 
Xdh(t): 



Xdh(t) = L + —!- — d 2 - \vriv uh (t 
2r]Vuh 2 \ 



■qvuh 



(19) 



Thus, according to Eq. ([Tq), R has a stationary shape 
(independent of time), but on a domain that extends 
downhill with time. Interestingly, we note that the mo- 
tion of x d h{t) is uniformly accelerated throughout the 
stage. This is a direct consequence of the non-linearity 
in Eq. ©. 



2. Central region 

In the central region, the boundary conditions have no 
influence on the evolution of R and h. The central re- 
gion is hence spatially defined by x u h(t) < x < x d h(t), 
and shrinks at both ends to disappear at the end of the 
stage. 

The evolution equation for R is the same as in the 
top region JEq. (|l7j)], but now we must impose the initial 
condition (|9j) (no boundary condition) . The solution (see 
Appendix) writes: 



R(x,t) = T)V uh t. 



(20) 



In the central region the rolling phase grows linearly 
with time and uniformly in space, thus forming a plateau 
(see Fig. |]-b). This constant growth rate is a conse- 
quence of the saturated form of the BCRE equations we 
have used [Eq. (g)]. The uniformity of the solution, on 
the other hand, stems from the fact that since none of 
the boundaries is at work, and since the initial static 
profile was uniform, the central region behaves like an 
infinite medium for which translational invariance is to 
be obeyed. 

We finally remark that solutions (|l8|) and (|2^) connect 
continuously at x — Xdh{t)- 
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3. Bottom region 

The bottom region is controlled by the bottom bound- 
ary condition, and spreads over the spatial interval 
< x < x u h (t) . In this region the evolution equation for 
R displays no more amplification (because x < x u h{t) = 
v uh t): dR/dt - T{R + d) dR/dx = 0. 

Since there is no constraint on R at the bottom of the 
cell, the condition on R is given by the physical assump- 
tion that it should be continuous across the border of the 
central and bottom regions, i.e. R(x — x u h(t), t) = r/v u ht 
for t > 0. 

This leads to the following expression for R: 



R(x, t) = - - (d+ - V v uht 



+ W\ (— + 1 - Tr > t ) +4— (* + r*). (2i) 

21 \ \V u h ) Vuh 

In this region also the height of rolling grains increases 
with time, due to an increasing input of material at the 
frontier with the central region. 



4- Derivation of time ti 

Stage I ends when the top and bottom regions meet, 
at t — t\ defined by x u h(t = t\) = Xdh(t = Using 
Eqs. ( |l4| ) and (|l9|), we easily obtain 



h = 




(22) 



As will be shown later, the maximum thickness of 
the avalanche, R m ax, is actually reached for t — t\ and 
x = x U h(t = ti) = Xdh{t = h). We can clearly see on 
Fig. [| that the i?-profile at time t = ti displays a cusp. 
As the prediction of the maximum amplitude R max is 
an i mport ant result of our analysis, we shall devote Sec- 
tion IV A to it, and defer the analytical derivation of 
Rmax and its application to physical examples until there. 

Fig. ||-b presents successive 'snapshots' of the rolling 
phase profile during Stage I. 



B. Stage II: The static profile reaches its final state 

Stage II starts at t = t\. At time t\, the two 'propa- 
gation fronts' of the boundary effects x u h(t) and Xdh{t) 
pass each other, and then pursue their respective mo- 
tions towards the opposite cell edge. Figure ^-a illus- 
trates this situation. As in Stage I, it appears that the 
cell is naturally divided in three spatial regions: a top re- 
gion [defined spatially (t) < x < L] , under the sole 
influence of the upper edge of the cell; a bottom region 
[0 < x < Xdh(t)], under the influence of the bottom edge; 



and finally, a central region [xdh(t) < x < x u h{t)], where 
in contrast with Stage I, the effects of both boundaries 
now combine. As another difference with the situation 
described in Stage I, the top and bottom regions progres- 
sively shrink, whereas the central one grows in extension 
(Fig. |a). 

Due to their motion, the fronts Xdh and x u h are bound 
to reach, sooner or later, the bottom and top end of the 
cell (respectively). At time ti [Eq. jlq)], the uphill front 
reaches the upper limit of the cell (x u h = L). The static 
profile is then in its relaxed final state, with a uniform 
slope 8 n . This is the end of Stage II, which is thus de- 
fined as the time interval t\ < t < t%. In most cases, as 
is discussed below, we expect the downhill front to reach 
the bottom edge before t = t%. 



1. Top region 

In this region, we have x > x u h, so that the evolu- 
tion equation ( |l7| ) still holds, and we still have to solve 
with respect to the upper boundary condition of Eq. 
Therefore, as in Stage I, R is given by Eq. (18) [but now, 
the lower limit of the domain on which this solution is 
valid is x u h{i), not Xdh(t)]. This top region shrinks, un- 
til finally disappearing when the uphill front reaches the 
upper end of the cell (t = £2)- 



2. Central region 

In this part, since x < x u h(t), there is no am- 
plification of the rolling amount, so that the right- 
hand side of the evolution equation of R vanishes: 
dR/dt - T(R + d) dR/dx = 0. Now, we further impose 
that R shall be continuous at the border with the top 
region, i.e.: 



R(x = x uh {t),t) = -d+ Jd 2 + 2(L - a:) 



r 



Solving these two equations together leads to look 
for R as one of the roots of the third-degree equa- 
tion R? + a 2 R 2 + ai(<)i? + a (x, t) = 0, with the follow- 
ing coefficients: 



«2 



( Vuh 

\ r 



3d 



ai(t) = 2^- [d + — r\(L - v uh t) 

a (x, t) = -2 V ^ {^(L - x) + 2d(L - v uh t)) . 
The solution, given by Cardano formulas [|22), writes: 



R( x ,t) = -^ + S-^, 



(23) 



with the auxiliary quantities [ 23 
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Q = ^(3di - a\) 
9a2ai — 27ao 



P 



2a 3 2 
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We saw in the previous section that the crest of 
the avalanche R m ax appeared at the end of Stage I. 
What happens to this crest during Stage II? It is easy 
to prove that the crest remains located on the down- 
hill front x c ih- Besides, x^h now moves at constant 
speed (in contrast with Stage I where it accelerated): 
Xdh{t) = v uh ti - T(R max + d)(t - h). 

We can also prove that the height of the crest remains 
constant (equal to R m ax) as it travels downhill, until it 
finally comes out of the cell. The exit time t ex it of this 
crest Rmax is obtained by solving Xdh(t) = 0: 



Vuhh 



r(i?„ 



d) 



1 



V u h 



T(R ri 



d) 



(24) 



Eq. (p3|). We present successive calculated 'snapshots' of 
the rolling amount R during this stage in Fig. &b. 



C. Stage III: The last grains are evacuated 

This stage lasts from t — t 2 until the end of the 
avalanche, at t = t en d- Both fronts have reached the 
edges of the cell at the end of Stage II {xdh = 0, x u h = L), 
and there is only one region (see Fig. More- 
over, as t > t 2l the slope of the static part is every- 
where n and no amplification of the rolling grains can 
take place; the rolling phase is simply convected down- 
wards. We now have to solve the evolution equation 
dR/dt — T{R + d) dR/dx = 0, with respect to the initial 
condition that Stage III evolves from what has been left 
by Stage II [i.e. R(x, t — t 2 ) as given by Eq. (23)]. 

Solving with the method of characteristics gives the 
following implicit solution 



R(x,t) = R(£,t 2 ), 



(25) 



where 



3. Bottom region 

The bottom region is defined as the region where 
< x < Xdhii)- The evolution equation for R is the same 
as in the central region, and we impose continuity at the 
border with the central region. We find that R(x,t) is 
given by Eq. (^l]) as in Stage I. Physically, in this region, 
we simply observe the convection of what was left in the 
bottom region at the end of Stage I. 

The bottom region disappears when the downhill front 
reaches the bottom end: Xdh(t) — 0, that is by defini- 
tion at time t — t ex it [Eq. (p4])]. To determine precisely 
the subsequent evolution of the avalanche, we must dis- 
cuss whether the disappearance of the bottom region oc- 
curs before the end of Stage II or not (that is, whether 
texit < h or t exlt > t 2 ). Using Eq. (g4|), we form the 
ratio: 

texit _ t\ A V u h \ 

~h~ ~ h. V T{R max + d)J ' 

Provided that the cell dimension L > 100<i and that 
V u h ~ Td (these requirements being usually satisfied 
for common experiments), we have TR max 3> v u h, and 
ti/t-z -C 1. Hence we generally expect t exit <C t 2 , and, 
consequently, as claimed earlier, in most cases the bot- 
tom region disappears before the end of Stage II. 

To resume, during Stage II, the central region extends 
both downhill and uphill with constant (though different) 
velocities at each end, progressively invading the whole 
cell. It reaches the bottom edge at time t ex u (in situ- 
ations where t ex u -C t 2 ), then the upper edge at time 
t 2 . At the end of Stage II, the central region occupies 
the entire cell, and R(x, t = t 2 ) is everywhere given by 



£ = * - T[R^,t 2 )+d](t-t 2 ) (26) 

(and < £ < L). 

The physical interpretation of these equations is ac- 
tually very simple: Eq. ( |25| ) states that the quantity of 
rolling species found in x at time t was previously located 
in £ at the beginning of Stage II {t — t 2 ). Equation ( pjjj ) 
gives a determination of this initial position £, by stating 
that from t 2 until the considered instant i, the quantity 
of grains moved with a constant speed r[i?(£, t 2 ) + d], de- 
pendent on the local height. In other words, during Stage 
III, the i?-profile left by Stage II is convected downhill, 
but each vertical slice rolls with its own velocity, which is 
a function of its height. The grains that were near the top 
edge of the cell at the end of Stage II are convected the 
most slowly, since there R was close to zero. The profile 
inherited from Stage II thus dilates upon rolling, under 
the effect of velocity inhomogeneities (Fig. |^-b) [ p4| . 

The last grains to fall off the cell are those that leave 
the top end of the cell at the beginning of Stage III, at 
time t 2 . Since at the top edge we have R = 0, these 
grains move with a constant speed vq = Td. At time i, 
they are located at xi ast (t) = L — Td (t — t 2 ), and the 
avalanche is extinct uphill: R — for x > xi as t (t) . 

Finally, the avalanche ends when the last grains reach 
the bottom limit of the cell (xi as t — 0), that is at time 
t end = t 2 + L/{Td). 

IV. DISCUSSION AND SIMPLE CHECKS 

A. Predictions for the maximum amplitude of the 
avalanche 



G 



1. Linear velocity profile 

Up to now, we have focused on flows displaying lin- 
ear veloci ty profiles. For such flows, as we saw in Sec- 
tion III A , the avalanche reaches its maximum amplitude 
Rmax at the end of Stage I, at time t = t\ [Eq. (p2|)]. 
The exact analytical expression of R m ax is easily found 
by using Eq. (|2(]) at time t\ (that is, the value of R given 
by the central region at the very moment it disappears) : 



Ri, 



Vuh 

r 



Vuh 

r 



2Liiv uh 

r ' 



(27) 



For large values of L, R max scales as: 



Rr; 



2 V V -fL, 



that is, as the square root of the system size L. 

Let us give a couple of numerical applications of this 
last expression. For the case of a standard laboratory 
experiment, with L = 1 m, d = 1 mm, v u h/T = 3d and 
77 ^ 0.1 rad, we find R ma x = 2.45 cm. In the case of a 
system at the scale of a desert dune, made of fine sand, 
we take L — 10 m, d — 0.2 mm and, with others param- 
eters unchanged, we get R m ax = 3.46 cm. One has to 
notice that R ma x is quite small, even for large systems 
as a sand dune. 

It is interesting to contrast this result with the work of 
Boutreux et al. [ |ll| , who carried the same calculation in 
an open cell configuration, but with a constant downhill 
convection velocity v(R) ~ v Q [instead of Eq. (Q)]. They 
found R m ax ~ "tlL. For the two above examples, this 
formula leads to maximum amplitudes of respectively 10 
cm and 1 m. The effect of the velocity gradient is thus to 
considerably limit the amplitude of avalanches, especially 
for large systems. 



2. Generalization for a power-law dependency 

In the beginning of this article, we quoted the work of 
Azanza et al. [ |l8| and of Pouliquen [ |l9| who find that 
the average speed of a chute flow of granular material on 
a rough plane is related to its thickness through a power- 
law relation v(R) ~ TR a with a close to 3/2. How- 
ever, as pointed out by Pouliquen [|lS[], the influence of 
the rough underlying bed plane on the rheology of chute 
flows is complex and not clearly understood, and might 
not be comparable to situations where the flow occurs 
on a free granular bed as has been considered in this pa- 
per. Since the question is still open, we will here present 
an intuitive derivation of R m ax valid for any power-law 
(undetermined exponent a) . To check the validity of this 
simple derivation, we first present it in the linear case 
a = 1, the generalization being then straightforward. 

Let us consider a point initially at the top edge of the 
cell. At t = 0, it starts being swept along by the granu- 
lar flow and we assume that this point travels with the 



local surface velocity of the flow v — TR. We are now 
interested in the temporal evolution of the rolling height 
R at this travelling point, which shall be computed from 
the Lagrangian derivative dR/dt = dR/dt + v dR/dx. 
As long as the amplification process takes place, we have 
with the use of Eq. (flfi|): dR/dt — 77 v u h- This implies 



R(t) = i]v uh t. 



(29) 



Hence, R(t) at the travelling point increases with time, 
as long as the amplification process lasts. After the am- 
plification has stopped, R at the travelling point keeps 
constant (since dR/dt = 0). Thus, R reaches its maxi- 
mum value Rmax at the end of the amplification. Let us 

(28) 

call t am p the duration of the amplification. We compute 



tn 



in the following way: the distance that the travel- 
ling point goes over during the amplification is of order 
~ L, so that tamp must verify: 



dx 



v dt 



i.e. L 



Tr) v uh tdt = -Trjv uh t am p 2 , 



in this calculation, we used v ~ TR(t)]. We finally find: 

gives the 



amp 



^2L/{Tnv uh ). 



Inserting this last expression into Eq. (|2 

Value Of Rmax- 



(30) 



This is exactly Eq. ( J28| ) found analytically, which 
was itself the limit of the complete expression of R m ax 
[Eq. (|2^)] for large values of L (greater than a hundred 
d), and v u h of order Td. 

The strongest assumption in the above simple 
derivation is that the amplification takes place 
over a distance ~ L. Rigorously, this distance is 
L — Xdh(t = ti); but what makes our simple derivation 
successful is that the position of the downhill front at 
time ti, Xdh(t = ti), is generally quite close to zero (for 
L greater than a hundred d and v u h of order Td) . 

We may now generalize the above results to a power- 
law dependency of the velocity v(R) ~ TR a . The same 
derivation leads us to the result: 



Rn 



(« 



r 



!/(«+!) 



(31) 



Note that R m ax diminishes as a increases. 

In particular, for a — 3/2, Eq. ( |3l| ) can be rewritten 
as Rmax ~ (5 nLv uh / 2 T) 2 / 5 . 
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B. Possible experimental checks 



2. Effects of polydispersity 



The loss of material at the bottom edge of the cell 
might be measured experimentally, and could be com- 
pared to the following theoretical prediction. This loss 
corresponds to the flow rate at the bottom of the cell 

r R(x=0,t) 



Q(x = 0,t)=tf 



v(z)dz, and is given by 



Q(x = 0,t) = -R(x = 0, i) 2 + Yd R(x = 0, t), (32) 



where R(x = 0, t) i s given by Eq. ( pl| ) during Stages I and 
II, and by Eq. ( pi5| ) during Stage III. Figure shows the 
predicted shape of Q(x = 0, t) as a function of time (solid 
curve). The curve displays a maximum at time t — t eX it, 
corresponding to the moment when the maximum ampli- 
tude R m ax rolls out of the cell. The maximum flow rate 
is obtained by replacing R(x 



0,t)byR max in Eq. fl32|). 
It is of interest to compare our prediction for the loss of 
material with that of Boutreux et al. [ [IlJ , who assumed 
a constant downhill velocity v in the rolling phase. This 
comparison, however, requires some caution: in our ap- 
proach, the granular flow is characterized by a constant 
velocity gradient T, whereas, in Boutreux et al., the de- 
scription is based on a typical d own hill convection veloc- 
ity of the grains v (see Section |l C|) . Figure compares 
the results of both approaches for the loss of material, 
assuming v ~ v u h — 3 Yd (see Ref. [ [Ti"|). 



C. Concluding remarks 



1. Regions of small R 



We notice that, during the avalanche, we had several 
spatial zones in the cell where R was close to (R < £'), 
e.g. at the upper edge of the cell, or at the end of the 
avalanche. Thus in these zones, the use of the saturated 
Eqs. (0) and (||) is not fully justified. In order to ob- 
tain a continuous description between the saturated case 
and the thin one, we could use the interpolated equations 
that have been proposed by de Gennes [ and studied 
in a model case by Boutreux and Raphael in Ref. [ : 



dh 
dt 

9R -Y(R 



RC dh 
R + ? dx 



d) dx- + ^R 



RC dh 



f dx 



The results of [ [r| show however that the physical 
behaviour is not dramatically changed, and that the de- 
scription in the zones of small R with saturated equations 
might be slightly wrong but qualitatively verified. 



It is of common knowledge that real granular materials 
are generally intrinsically polydisperse. This may have 
drastic effects on the behavior of the flow, and capturing 
more precisely the physics of real avalanches would cer- 
tainly suppose to take polydispersity into account. How- 
ever, the treatment of full polydispersity is a difficult 
task. Yet, the BCRE equations have been extended to 
the case of binary mixtures [ |2|,|2(| , and it could be in- 
teresting to study the changes brought up in this case by 
a velocity gradient in the flow. 



3. Domain of validity of the BCRE approach 

The general approach introduced by Bouchaud et al. 
to describe surface flows is rather phenomcnological, and 
as pointed out by Bouchaud and Cates [ [l4| , we still 
lack criteria to determine the range of physical situa- 
tions to which it can be successfully applied. In a recent 
work, Douady et al. [ [l6| proposed a justification of the 
BCRE modelization on the basis of hydrodynamic con- 
servation laws. According to these authors, the form of 
the BCRE equations should not remain invariant when 
different laws are chosen for the velocity profile in the 
flow. Douady et al. argue that only for a velocity linear 
in R (or constant) shall the equations take the simple 
form of our equations [Eqs. (^) and (||)]; in other cases, 
they find that a supplementary term coupling R and h 
should add in Eq. (|7|). Certainly, more work needs to 
be done in this direction in order to exactly assert the 
domain of validity of the BCRE analysis. 
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APPENDIX: 
1. Method for solving the evolution equation of R 

Eq. (|l6|) is a first-order partial differential equation, of 
the quasi-linear class, that is, linear in the first deriva- 
tives. Such equations can be solved by the well-known 
method of characteristics. See, for example, Ref. [ |2lf|. 

More specifically, we will solve Eq. ([l6]) along 
characteristic curves given in the parametric form 
{t(s),x(s), R(s)}, with s the parameter. The functions 
i(s), x(s) and R(s) are derived from the set of coupled 
ordinary differential equations: 
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dt _ ^ 
ds 

^=-T(R + d) 
ds 

— = T]V uh H{x - v uh t). 

ds 



(Al) 



By integration, one founds the equations for the char- 
acteristics with unspecified integration constants. One 
then imposes the boundary and/or initial conditions to 
identify these constants. 

We here give the detailed calculations only for the first 
stage of the avalanche. The derivations are separated 
into the different spatial regions that were defined ear- 
lier, and we will show how they naturally emerge from 
the derivations. 



2. Top region 



In this region, Eqs. (Al) become: 



dt 

ds 
dx 

ds 
dR 

ds 



-T(R+d) 



(A2) 



V v uh- 



We also use the boundary condition R(x — L,t) = 
for t > 0, which we parametrize with the parameter £. 
For simplicity's sake, on each characteristic crossing the 
boundary curve we arbitrarily choose the value of the 
parameter s to be zero at the crossing point. This deter- 
mines the integration constants to be: 



t(s = 0) = £ 
x(s =Q)=L 
R(s = 0) = 0. 



(A3) 



Note that £ > 0, since t > (the experiment started 
at time t = on). 



Solving for Eqs. ( |A2| ) together with Eqs. ( |A3| ) gives the 
equations of the characteristic curves: 



t(s) 



x(s) 
R(s) = rjv uh s 



(A4) 
(A5) 
(A6) 



We now want to write the solution R explicitly in t erm s 
of x and t, so that we have to eliminate £ and s. Eq. (A5) 
can be so lved to give s as a function of x, and replacing 
into Eq. (A6) brings the analytical solution: 



R(x,t) 



d 2 + 2(L-x 



rjv uh 



(A7) 



We now have to verify th e co ndition that £ > 0. By 
combining Eq. (A4) with (A6) this condition can be 
rewritten as t > R{x^t)/r]v u h. Replacing R into this 



inequal ity [by ( A7)j gives us a spatial condition for so- 
lution (A7) to be valid: we must have x > Xdh{i), where 



x d h{t) = L + 



1 



Fr)v uh (t ■ 



2vvuh 2 r/v uh 

This is the mathematical origin of the 'downhill front' 
that we described intuitively in the main text as the limit 
of extension of the boundary effects originating in the up- 
per edge of the cell. 



3. Central region 

The evolution equation for R is the same as in the top 
zone, so that the differential equ ation s giving the charac- 
teristics also are the same [Eqs. (A2)]. But now we must 
impose the initial condition R(x,t = 0) = 0, which gives 
the following set of initial conditions for the characteris- 
tics: t(s = 0) = 0, x(s = 0) = £, R(s = 0) = (and 
< £ < L). We obtain: 



t(s) - 8 

x(s) = --Triv uh s 2 - Yds 
R(s) = rjv uh s. 



L 



(A8) 
(A9) 
(A10) 

Combining (A8) and ( A10 ) gives an explicit solution 
for R: R{x, t) — T]v u ht. 

This solution is valid in a certain spatial domain. It 
is limited upwards by the top region [i.e. x < Xdhit)]. 
It is also limited downwards by x u h(t), because at this 
point the form of the evolution equation of R changes 
(the amplification term vanishes) , and consequently does 
the form of the differential equations that give the char- 
acteristics. 



4. Bottom region 



In this region, Eqs. (Al) are given by: dt/ds = 
1, dx/ds = —T(R + d), dR/ds = 0. 

Here, the boundary condition is given by the continuity 
of R at the border of the central and the bottom zones: 
R(x — x u h(t), t) = rjVuht for t > 0. This gives the initial 
conditions: t(s = 0) = £, x(s = 0) = v u h£, R(s = Q) = 
V v uh£,- Solving and rewriting R explicitly in terms of x 
and t leads to the solution: 

Vuh 

2 v ' r 



R{x,t) 




- rjv uh t 



l-T-qt 



4—(x + Tdt), 

Vuh 



valid for < x < x u h(t). 



which is the same as Eq. (p 
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FIG. 1. The basic assumption of the BCRE picture is that there is a sharp distinction between immobile grains with a profile 
h(x,t), and rolling particles with a local amount R(x,t). The immobile grains constitute the 'static phase' and the rolling ones 
the 'rolling phase'. The local slope of the static profile is called 8(x,t). 

FIG. 2. Example of an open cell, so as to let the rolling material flow out. We suppose that the avalanche starts precisely at 
P — Umax (see text). 

FIG. 3. The profile of the static grains for < t < t 2 . At the left of x uh (t), the slope has relaxed to its final value 9 n . On 
the right, it still has the initial angle max . 

FIG. 4. (a) Position (dotted lines) and motion (arrows) of the 'downhill' and 'uphill' fronts during Stage I. The respective 
sizes of the static phase (dark) and rolling phase (light) have been modified for clarity purposes. The positions of the fronts 
naturally define three regions, with specific physical meaning (see text): bottom (I), central (2), and top (3). (b) Evolution of 
the rolling phase in the cell during Stage I. The plot presents R vs. the position x, at successive dates (R and x are given in 
grain diameters d, and the parameters are v u h = STd, 77 = 0.1 rad and L — 1000(1 Note the different horizontal and vertical 
scales). The amount of rolling grains grows with time in the whole cell. Regions borders correspond to slope discontinuities. 
For t = ti, the profile presents a cusp, where the maximum thickness R ma x of the avalanche is reached. 

FIG. 5. (a) 'Downhill' and 'uphill' fronts during Stage II. The fronts naturally define three spatial regions: bottom (1), 
central (2), and top (3). (b) Evolution of the rolling phase in the cell during Stage II. The plot presents R vs. the position 
x at successive dates from t = t exit to t = t 2 (R and x are given in grain diameters d; v uh , r\ and L as in previous figures. 
Note the different horizontal and vertical scales). The amount of rolling grains globally decreases in the cell. Regions borders 
correspond to slope discontinuities. At t = t ex it, the peak amplitude R ma x reaches the bottom end of the cell. 

FIG. 6. (a) Both fronts have reached the cell borders (in most cases; see text), there is only one region in the cell, (b) 
Evolution of the rolling phase in the cell during Stage III (R and x are given in grain diameters d; v u h, i) and L as in previous 
figures). The profile at the beginning of the stage (t = t 2 ) is progressively convected downhill, but dilates at the same time, 
as seen on the plots at t = £2 + 0.05(t en d — t 2 ) and t = t 2 + 0.2(t en d — t 2 ). This is because thicker vertical slices roll faster 
than thinner ones. We see that the grains at the top edge of the cell (x = L) at time t = t 2 are the last ones to evacuate; the 
avalanche is extinct uphill (R — 0). 

FIG. 7. Loss of material at the bottom edge of the cell as a function of time. The loss is given by the flow rate Q(x — 0,t) 
(Q is in units of Td 2 , t in units of ~ (d/g) 1 ^ 2 ; v u h, r\ and L as in previous figures). Solid line: predicted shape with a linear 
velocity profile in the flow; dashed line: predicted shape in the case of a constant velocity profile in the flow, from Boutreux et 
al, with the choice v = v uh (see text). 
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